A STOCHASTIC MODEL FOR BACTERIOPHAGE THERAPIES 
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Abstract. In this article, we analyze a system modeling bacteriophage treatments for 
infections in a noisy context. In the small noise regime, we show that after a reasonable 
amount of time the system is close to a sane equilibrium (which is a relevant biologic 
information) with high probability. Mathematically speaking, our study hinges on con- 
centration techniques for delayed stochastic differential equations. 



o 

CN . 1. Introduction 

^| In the last years Bacteriophage therapies are attracting the attention of several sci- 

C/} ■ entitle studies. They can be a new and powerful tool to treat bacterial infections or to 

\Q . prevent them applying the treatment to animals such as poultry or swine. Very roughly 

speaking, they consist in inoculating a (benign) virus in order to kill the bacteria known 
fV| to be responsible of a certain disease. This kind of treatment is known since the beginning 

Oh ■ of the 20 th century, but has been in disuse in the Western world, erased by antibiotic ther- 

apies. However, a small activity in this domain has survived in the USSR, and it is now 
re-emerging (at least at an experimental level). Among the reasons of this re-emersion 
we can find the progressive slowdown in antibiotic efficiency (antibiotic resistance). Re- 
ported recent experiments include animal diseases like hemorrhagic septicemia in cattle or 
atrophic rhinitis in swine, and a need for suitable mathematical models is now expressed 
^ | by the community. 

Let us be a little more specific about the (lytic) bacteriophage mechanism: after attach- 
ment, the virus' genetic material penetrates into the bacteria and use the host's replication 
mechanism to self- replicate. Once this is done, the bacteria is completely spoiled while 
^ . new viruses are released, ready to attack other bacteria. It should be noticed at this 

point that among the advantages expected from the therapy is the fact that it focuses 
on one specific bacteria, while antibiotics also attack autochthonous microbiota. Roughly 
speaking, it is also believed that viruses are likely to adapt themselves to mutations of 
their host bacteria. 

At a mathematical level, whenever the mobility of the different biological actors is 
high enough, bacteriophage systems can be modeled by a kind of predator-prey equation. 
Namely, set St (resp. Q t ) for the bacteria (resp. bacteriophages) concentration at time t. 
Consider a truncated identity function o : R + — > IR + , such that a G C°°, cr(x) = x 
whenever < x < M and a(x) = M + 1 for x > M + 1. Then a model for the evolution 
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of the couple (S, Q) is as follows: 

dS t = [a - k<j(Q t )\ S t dt 

dQ t =[d- mQt - ka(Q t )St + kbe^a(Q t ^)S t ^} dt, 

where a is the reproducing rate of the bacteria and k is the adsorption rate. In equation 
P]), d also stands for the quantity of bacteriophages inoculated per unit of time, m is their 
death rate, we denote by b the number of bacteriophages which is released after replication 
within the bacteria cell, ( is the delay necessary to the reproduction of bacteriophages 
(called latency time) and the coefficient represents an attenuation in the release 

of bacteriophages (given by the expected number of bacteria cell's deaths during the 
latency time, where /i is the bacteria's death rate). A given initial condition (Sq,Qq) is 
also specified. When modeling biological phenomena, one usually assumes models like (CQ), 
where a is replaced by the identity function. We have considered here the truncation of the 
identity a in order to manipulate bounded coefficients in our equations, but our parameter 
M can also be interpreted as a maximal infection rate of bacteria by bacteriophages. 
One should also be aware of the fact that the latency time ( (which can be seen as the 
reproduction time of the bacteriophages within the bacteria) cannot be neglected, and is 
generally of the same order (about 20mn) as the experiment length (about 60mn). 

According to the values of the different parameters of the system and of the initial 
conditions, different types of equilibriums for equation (CQ) might emerge. We shall focus 
in the sequel on the simplest of these regimes, namely when d is large enough (the exact 
condition is kd/m > a). This makes the mathematical analysis easier, and it corresponds 
to the existence of a unique stable steady state Eq = (0, d/m) for our system (in particular 
bacteria have been eradicated). Notice however that we can perfectly assume the regime 
kd/m > a since the treatment allows to inject high quantities of viruses. One should 
also mention a natural generalization of our problem: Consider the action of several 
varieties bacteriophages, which is an option widely considered among practitioners. We 
have restricted our analysis here to a simplified situation for sake of readability. 

It is perfectly assumable that noise will appear when collecting data from laboratory 
tests. Moreover, when one wishes to go from in vitro to in vivo modeling, it is commonly 
accepted that noisy versions of the differential systems at stake have to be considered. 
This program has been carried out e.g. for HIV dynamics in |3] and for bacteriophages 
in marine organisms in [2]. In those references it is always assumed that the noise enters 
in a bilinear way, which is quite natural in this situation and ensures positivity of the 
solution. We shall take up this strategy here, and consider system ([I]) with a small 
random perturbation of the form 

dSf = [a- ka{Ql)\ Stdt + ea(St) o dW? 

dQ\ =[d- mQt - ka(Q e t )St + kbe-^a{Q £ t _ c )S £ t _ c ] dt + ea(Q e t ) o dWf , 

where e is a small positive coefficient and W = (W 1 ,W 2 ) is a 2-dimensional Brownian 
motion defined on a complete probability space (Q, J 7 , P). Our aim will then be to prove 
that for a time To within a reasonable range, the couple Z e TQ '■= (S £ Q ,Q^. ) is not too far 
away from its stable equilibrium Eq. Note that reasonable range is meant here as a time 
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dQ £ t = [d- mQ £ t + k(b - l)a(Ql)S E t ] dt + ea(Q £ t ) o dW? 



which corresponds to the order of both the latency delay and the time when the immune 
system of the animal can cope with the remaining bacteria. 

As we shall see in the sequel, the treatment of equation (j2J) involves the introduction 
of some rather technical assumptions on our coefficients. For sake of readability, we have 
thus decided to handle first the following system without delay: 

dSt = [a- ka{Ql)] Stdt + ea(St) o dW} 

t > 

where we notice that the only difference between ([2]) and ([3]) is that we have set £ = in 
the latter. 

The main advantage of equation fl3]) lies into the fact that we are able to work under 
the following rather simple set of assumptions: 

Hypothesis 1.1. We will suppose that the coefficients of equation (TJ)] satisfy: 
(i) The initial condition (So, Qo) of the system lies into the region 

mM - d 



0. 



x [d/m } M\. 



k(b-l)M_ 

(ii) The coefficient 7 = kd/m — a is strictly positive and M > d/m. 

We shall also use extensively the following notations: 

Notation 1.2. The letters c, ci,c%,... will stand for universal constants, whose exact 
value is irrelevant. For a continuous function f , we set ||/||oo,/ = sup^gj 

Then the previous loose considerations about convergence to Eq can be summarized in 
the following theorem, which is the main result of our paper for our bacteriophage system 
without delay: 

Theorem 1.3. Given positive initial conditions, equation (TJ|) admits a unique solution 
which is almost surely an element o/C(R+, Ri). Assume furthermore Hypothesis set 
77 = m/2 A 7 and consider 3 constants 1 < K\ < K2 < ^3- Then there exists po such that 
for any p < po and any interval of time of the form I = [ki ln(c/p)/f], K2 hi(c/p)/7?], we 
have 

P(\\Z S - EoW^j > 2p) < exp (-^-) , ( 4 ) 
where X is a constant satisfying A > K3/77. 

Remark 1.4. Relation (j4j) can be interpreted in the following manner: assume that we 
observe a noise with intensity e. Then the kind of deviation we might expect from the 
noisy system ([3]) with respect to the equilibrium Eq is of order with $ = 277/^3. This 
range of deviation happens at a time scale of order ln(p _1 )/?7. 

A second part of our analysis is then devoted to the more realistic delayed system, for 
which we obtain a result which is analogous to Theorem 11.31 

Theorem 1.5. Equation fl4]) still holds for the delayed system under some slightly 
more restrictive conditions on the initial condition which shall be specified at Hypothe- 
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Theorem 11.51 can be seen as the main result of the current paper, and deserves some 
additional comments: 

(1) We have produced a concentration type result instead of a large deviation principle 
for equation ([3]), because it seemed more adapted to our biological context. Indeed, in the 
current situation one wishes to know how far we might be from the desired equilibrium at 
a given fixed time, instead of producing asymptotic results as in the large deviation theory. 
At a technical level however, we rely on large deviation type tools, and in particular on 
an extensive use of exponential inequalities for martingales. 

(2) Let us compare our result with [21 E], which deal with closely related systems. The 
interesting article [2] is concerned with a predator-prey system similar to ours, but it 
assumes that a linearization procedure around equilibrium in the highly nonlinear situa- 
tion ([3]) can be performed. The analysis relies then heavily on this unjustified step. As 
far as [3] is concerned, it roughly shows that if the noise intensity of the system is high 
enough, then HIV epidemics can be kept under control (in terms of exponential stability). 
This is a valuable information, but far away from our point of view which assumes a 
low intensity for the noise. We should also mention the related thorough deterministic 
studies pQ El E] 

(3) Mathematically speaking, it would certainly be interesting to play with the rich 
picture produced by equation ([I]) and its perturbed version in terms of stable and instable 
equilibrium. We have not delved deeper into this direction because it did not seem directly 
relevant to the biological problem at stake. It should be pointed out however that the 
analysis of our random dynamical system (EJ) is non standard due to the coefficient d, 
which accounts for the bacteriophage inoculation. Many of our considerations below will 
be devoted to handle this problem. 

Our article is structured as follows: Section [2] is devoted to some preliminary con- 
siderations (existence and uniqueness results for our stochastic systems, convergence to 
equilibrium for the corresponding deterministic equations). Then we show our concen- 
tration results at Section El Finally some simulations are lead at Section 2] in order to 
illustrate the theoretical results. 

2. Preliminaries 

In this section, we give some basic results concerning our competition system. We first 
establish existence and uniqueness for the solution to the perturbed system ([2]), starting 
from the simpler system ([3]). Then we deduce some properties for the equilibria of the 
deterministic counterpart of both systems ([2]) and ([3]). 

2.1. Existence, uniqueness and positivity of solution. Recall that we are consid- 
ering the perturbed problem ([2]), with a coefficient a and some initial conditions of the 
following form: 

Hypothesis 2.1. The coefficients of our differential systems satisfy the following assump- 
tions: 

(i) The function a : M + — > R + is such that o G C°° , and satisfies a(x) = x for < x < M 
and a(x) = M + 1 for x > M + 1. We also assume that < cr'(x) < C for all x G M +7 
with a constant C such that C > 1. 
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(ii) As far as the initial condition is concerned, we assume that it is given as continuous 
positive functions {So, t ,Qo,t', —( < t < 0}. In case of the non delayed system (pQ) ; it is 
simply given by two positive constants [Sq,Qo). 

Due to the fact that we have assumed a bounded coefficient a, the existence and 
uniqueness of the solution to our differential system is a matter of standard considerations. 

Theorem 2.2 (Global existence of solution). For any positive initial condition there exists 
a unique solution of §Z§, which is defined for all t > 0. 

Proof. It is readily checked that the coefficients of the equation are locally Lipschitz with 
linear growth. The existence and uniqueness of the solution is then a direct consequence 
of classical results (see e.g. [6J Section 5.2] for the non delayed system and [H] for the 
delayed one). 

□ 

Positivity of the solution is also an important feature, if we want the quantities St, Qt 
to be biologically meaningful. Moreover, part of our analysis will rely on this property, 
that we label for further use: 

Proposition 2.3 (Positivity). If we take positive initial conditions Sqj > 0, Qo,t > for 
all t G [— C, 0] for the system (j2J), then the solution fulfills Sf > 0, Q\ > for allt > 0. 

Proof. Let us first consider the system with o~{x) = x for all x, namely: 

dSt = [a- kQt] S?dt + eSt ° dW, 1 

dQl =[d- mQl - kQtSt + k b e'^Q^S^] dt + eQ\ o dWf, 

with initial condition {So tt ,Qo,t)- Assuming existence and uniqueness of the solution to 
(jSJ), we shall prove that S|, Q\ > for all t > almost surely. 

Indeed, after the change of variables x t = e~ EW * Sf, yt = e~ eW tQf, we can recast ([5]) 
into the following system of differential equations with random coefficients: 



x' t = (a - ke £W ?y t ^ x t 



(6) 

y > = de~ eW ? - my t - ke eW *x t y t + k be^-^^-^-Oyt^xt-c 

with initial conditions x°(t) = Sqj > 0, y°(t) = <3o,t > for all t G [—£,()]. Then, the 
positivity of x t is immediate from the representation 

ft 



x t = x°(0) exp | jf (a - ke eW °y s )ds^ > 0. 



In order to see the positivity of y t let us observe that for y°(0) = we have y' = 
d+ A;6e" MC_e(H/ ° 2 " H/ -c" w -c ) ?/_ ? x_ ( > 0. Therefore, for all initial condition y > there 
exists 5 > such that y t > for all t G (0,5). Let us suppose now that y t < for 
some t > 0, and let to = inf{t > | yt < 0}. Due to the continuity of the solution we 
have that y to = 0. Then y' to = de~ eWt o + kbe~^~ e ^ Wt o~ Wt o-c~ Wt o-c' > y to _^Xt -( > 0, which 
is impossible since it would yield y t > for t G [to, to + 5) for 5 small enough. This 
contradiction means exactly that y t > for all t > 0. 



6 



X. BARDINA, D. BASCOMPTE, C. ROVIRA, AND S. TINDEL 



Now that we have the positivity for system (jSJ), we can prove the positivity for ([2]) in 
the following way. Let us first handle the case of 5f , and assume that the initial condition 
is such that 5o,o > M. Set then rj^ 5 = inf{t > such that 5f < M/2}, and observe that 
s 1S a stopping time for the natural filtration of the Brownian motion W, such that S e 
has remained positive until t® 4 s . Furthermore, the strong Markov property for (S £ ,Q e ) 
entails that the process 



also satisfies (j2J) on the set Qm,s = G f2; Tj^s < oo}, with an initial condition 5*0,0 = 
M/2. With these considerations in mind, we can assume that the initial condition of our 
differential system satisfies 5*0,0 < M. 

With such an initial condition we can conclude the positivity of Sf until the stopping 
time Tms = inf{t > such that Sf > M} as we have done for the system (JSJ), since up 
to time Tms we have cr(5f) = 5f. Then, invoking again the strong Markov property, 
we can also guarantee positivity until time t]^ s = inf{t > Tm s sucn that Sf < M/2} as 
above. We are now in a position to obtain the positivity of Sf until time f\ 4 S = inf{£ > 
t m s sucn that Sf > M}, once again with the same reasoning than for the system (JSJ). 
The global positivity of Sf on any interval of the form [t^ s , T^/g] for k > now follows 
by iteration of this reasoning. 

It remains to show that lim^oo s = oo. This is easily obtained by combining the 
following two ingredients: 

(i) The increments {r^/g — s ; k > 0} form a i.i.d sequence by a simple application of 
the strong Markov property. 

(ii) Owing to the specific coefficients we have for equation (j2J), it can be checked that for 
any 772 > one can find 771 > small enough such that ~P{t\j s > 771) > 1 — 772. Details of 
this assertion are omitted for sake of conciseness. 

We let the reader check that the positivity of Q\ can be obtained along the same lines, 
which ends the proof. 



Remark 2.4. Using the a priori positivity properties stated above, we could have also 
obtained existence and uniqueness of the solution for system (j5J). We did not include 
those developments for sake of conciseness. 

2.2. Analysis of the deterministic non delayed system. This section is devoted to 
the analysis of the non perturbed system corresponding to . Namely, we shall consider 
the following dynamical system: 



We will give some sufficient conditions for the existence of a unique stable equilibrium E$ 
and then show exponential convergence to this equilibrium. 

Let us start with the basic results we shall need about equilibria of (J2J). 




□ 



dS t = [a — ka{Qtj\ S t dt 

dQ t = [d- mQ t + k(b - l)cr(Qt)S t ] dt. 



(7) 
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Theorem 2.5. If either M + 1 < |f or M > § and — > a, system (J7|) /ias a unique (posi- 
tive) steady state E = (0, -^). Moreover, the bacteria-free equilibrium E is asymptotically 
stable for — > a and M > — . 

Proof. To obtain the equilibria, we have to find the solutions of the following equation: 

= (a-ka(Q))S 

, (8) 
= d-mQ + k(b-l)a(Q)S, 

where S, Q are positive constants. 

Owing to the first equation we have either S = or a — ka(Q) = 0. Since S = 
and the second equation imply Q = — , we have that bacteria-free equilibrium E exists 
for any value of the parameters. In the case M + 1 < f one can observe that no other 

equilibrium exists (since a — ka(Q) > for any Q). 

Taking M > f , a — ka(Q) = if and only if Q = |f . Then, using the second equation 
in (jSj), we have 

t a /, i\ A a ma — kd 

Q = d-m- + (b-l)aS =► 5 = ^-^, 

which is positive only for a > So we have proved the first part of the result. 
For the second part, the Jacobian matrix of system (J2J) at i?o is 

a-Jfetr^) 
v fc(6-l)a(i) -m y 

The eigenvalues of this matrix are easily shown to be Ao = a — ka(-) and Ai = — m, 
which are negative for — > a and M > — . 

□ 

We now wish to study the rate of convergence towards the E Q equilibrium in the stable 
case (i.e., when kd/m > a and M > —). The main result we obtain to this respect is: 

Theorem 2.6. Under Hypothesis li.il and \2.1\ the solution of system (j7|) with initial 
condition 

mM - d 

0, 



A :-- 



x [d/m, M] 



(S ,Q )e rk{b _ 1)M 

exponentially converges to the equilibrium E : 

\(S t ,Q t )-E \<ce^\ Wl th r] = j A — , (9) 

where we recall that 7 = — — a > 0. 

Proof. In order to prove our claim, we first have to show that the region R := [0, j^zjjfj] X 

[^,M] C [0,M] 2 is left invariant by equation (0). Towards this aim, we can invoke the 
same method we will use in Proposition 12. S\ and we let the reader check the details. 

Now, since we have Q t < M for all t, we can consider a(x) = x in equation ([7j). We 
will consider a version of this system centered at E by means of the change of variables 
S = S,Q = Q — d/m. This leads to the system 
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SI 



kQtSt 
kd 



(10) 



-J s t 

Q' = - m Q t + -(b _ 1)^ + k (b _ l)Q tSt . 

m 

Notice that, according to our set of assumptions concerning the initial conditions, we have 
Sq > and Qo > 0. Thus the solution to (|10|) will remain positive for all t > (it can be 
deduced from R being invariant, or can be proved just like in Proposition 12.3 j) . 

Now, from the first equation in (1101) . we have that S' t < —jSt. This implies S t < Soe -7 *, 
proving that St exponentially converges to zero. 

Owing to the second equation in ( flT)]) and using positivity properties of the solution, 
we also get 



Q't < -mQt + k(b - l)S e" 



-it 



d 
m 



Q, 



Finally, the variation of constants method will lead to the stated result, following the 
same steps we will detail later in the proof of Theorem 12.101 

□ 

2.3. Analysis of the deterministic delayed system. We now try to generalize the 
results of Section |2~21 to our deterministic delayed system (JTJ). To this aim, we shall work 
under the following assumptions. 

Hypothesis 2.7. We will suppose that the coefficients of equation satisfy the following 
conditions, valid for any t G [— 0]: 

(i) The initial condition (So,*, Qo,t) °f the system lies into the region 



R := [0, M] x 



-,M 
m 



;ii) We have be-^Q 0>t S 0it > ±S , , and fee"* > 1. 
^iii) The condition Sqj < k T A1 ^c'L is satisfied. 



A first step towards exponential stability is then the invariance of a certain region under 
our dynamical system: 

Proposition 2.8. Under Hypothesis \2.1\ and \2.7[ the region 



R :-- 



0. 



mM -d 
kbe-^M 



d 



m 



M 



C [0,M] 2 



is left invariant by equation ([T]). 

Proof. We separate the analysis of S and Q in two steps. 

Step 1: boundedness of S. Since S is obviously positive (along the same lines as for 
equation (jSJ)) and owing to the fact that S' t = (a — ka(Q t )) S t we obtain that 



St < whenever Q t > — , and St > whenever Q t < — . 

k k 



Furthermore, our system starts from an initial condition Qo,o ^ 
increasing as long as Q remains in the interval [— , oo). 



> 



Thus S is non 
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Let us now observe what happens in the limiting case Qo,o = ~ : recalling that our 
initial conditions are denoted by So,*, Qo,t for t G [— £, 0], we have 

Q' = -k-S 0fi + kbe-^a(Qo,- ( )S 0t . c = k (V^Q ,-c*Vc - -S 0>0 ) > 0, 

where we have used the fact that &e _/ ^(5o,-£>S'o,-c > m^ofl- According to this inequality, 
we obtain the existence of a strictly positive e such that Q t > — for all t G (0, e). We 
thus introduce the quantity to = inf{t > : Qt — — }, and notice that we have 

Q' to = -k±S t0 + kbe-^a(Q t0 _ ( )S t0 _ c . 

We can now distinguish two cases: 

(1) If to > C? since St is non-increasing in [0,t ], ^to-C — ^to an d hence 

, d 



Q' to > kS t0 (be-"a(Q 



t0 - a I > 0, 

m 



due to the fact that be~**> > 1, M > — and Q tn -c > — • 
(2) If t < C> since S to < S 0)0 we obtain 

Q't > -k^S 0)0 + fcfoe _K (r((5o ! to-c)' S 'o ! to-C = fc f 6e_K( 5o,to-C<So,t -c ~ ""^J > °' 

where we have used the fact that be~^Qo ; tSo,t > 5*0,0 f° r all t G [— C, 0]. 
This discussion allows thus to conclude that to cannot be a finite time. Indeed, we should 
have Q[ > and hence Q increasing in a neighborhood of to, while Q should be decreasing 
in a neighborhood of t according to its very definition. We have thus reached the following 
partial conclusion: 

Qt > — , t h> St decreasing, St > 0. 
m 

In particular, any interval of the form [0, L] for L > is left invariant by t y St- 

Step 2: boundedness of Q. Our claim is now reduced to prove that for (So t t, Qo,t) G R we 
have Q t < M for all t > 0. 

To this aim notice that, whenever Qo,o = M we have 

Q' Q = d-mM -kMSofi + kbe-^aiQo^So^ 

< d-mM + kbe' K MS - C < 0, 

where we recall that So,-c < according to Hypothesis 12. 71 This yields the existence 

of e > such that Q t < M for all t G (0, e). 

We now define t\ = inf {t > : = M}. It is readily checked that 
Q' h = d-mM -kMSt.+kbe-^o-iQt^St.-c 
= d - mM - kMS tl + kbe^Qt^St.-c 

< d-mM + kbe~ K MS tl - 
and we can distinguish again two cases: 
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(1) If ti > C, thanks to the fact that 1 H y S t is non-increasing on [0,ti], we have 

Q' tl < d - mM + kbe~ K MS 0fi < 0, 

since we have assumed that S^o < ir^^Tj ■ 

(2) If tt < ( then 

Q' tl <d- mM + kbe-^MS 0ttl -c < 0, 

thanks to the fact that S 0>t < t f^, for all t G [-£, 0]. 

As for the discussion of the previous step, this allows thus to conclude that t\ cannot be 
a finite time, due to the contradiction Q' t < and Q t < Q tl for all t G (0, t\). We have 
thus shown Q t < M for all t > 0, which finishes the proof. 

□ 

Remark 2.9. Before stating the exponential convergence to the bacteria-free equilibrium 
result, let us observe that Theorem 12.51 still holds true for the delayed system ([TJ). It can 
be easily checked using exactly the same steps we have done for the non-delayed system. 

We are now ready to state our result on exponential convergence of the delayed dynam- 
ics: 

Theorem 2.10. Assume Hypothesis \2.1\ and \2. 71 are satisfied, and let R be the 

region defined at Proposition \2.£\ Then the solution of system (CQ) with initial condition 
(So,Qo) G R exponentially converges to the equilibrium E : 

in 

\{S t ,Q t )-E \ <ce-*, with 77 = 7 A-, (11) 
where we recall that 7 = — — a > 0. 

Proof. According to Proposition 12.81 we have Q(t) < M for all — ( < t < 00 under our 
standing assumptions. Hence one can recast equation (TJQ) as 

J dS t = (a — kQ t ) S t dt 

\ dQ t = (d- mQ t - kQ t S t + kbe-rtQt-cSt-c) dt 

Let us perform now the change of variables Q = Q — — . This transforms the previous 
system into 

dS t = ( a + S t dt 

dQ t = (d - m{Q t + i) - k{Q t + + kbe-^(Q t - C + £)S t -^ dt. 
Equivalently, our new system is: 

dS t = - (jS t + kQ t S t ^) dt 

dQ t = (-mQ t - k±S t - kQtSt + k£be~^S t _ c + kbe^Q t ^S t ^) dt. 

Observe now that Proposition 12.81 asserts that Q t > — for all t > 0, which means that 
Qt > 0. With our change of variables, we have also shifted our equilibrium to the point 
(0, 0). We now wish to prove that St and Qt exponentially converge to 0. 
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The bound on St is easily obtained: just note that 

dS t < -7<S*dt, 

which yields St < S^o e -7 *. As far as Q t is concerned, one gets the bound 



7 (t-C) 



^1 < _ m g t + k—be-i* 5 ,o e" 7 ^ + kbe^Q t ^S „ e" 
dr m 

< -mQ t + fc6e- K ^o o e" 7 ^ f — + M - -] 

\m m J 

= —mQ t + ce~ 7 *, 

with c = kbMSo : o e^ 7- ^, and where we have used the fact that Q t < M uniformly in t. 

Invoking now the variation of constant method, it is readily checked that equation 
it = —mxt + ce" 7 * with initial condition xq = Qo,o can be explicitly solved as 

x{t) = e- mt ( Q ,o + — ^— (e (m ~ 7)t - 1) 
\ m — 7 

Qo,o ~ -?—) e~ mt + e" 7 *. 
m — 7/ m — 7 

By comparison, this entails the inequality < Ci e _,?t , where Ci = max(C)o,o — T^ji 
and 77 = m A 7. Our proof is now finished. 

□ 

3. Fluctuations of the random system 

Let us summarize the information we have obtained up to now in the non delayed case: 
we are considering the system 

dSt = [a- ka(Qt)} Stdt + ea(St) o dW t ' 

dQ £ t = [d- mQ e t + k{b - l)a{Ql)Sf[ dt + ea{Q e t ) o dW?. 

Under Hypothesis 11.11 and 12. 1\ we have shown the existence of a unique equilibrium E = 
(0, d/m) for the deterministic system (J7j), corresponding to (|12|) with £ = 0. Furthermore, 
we have constructed a region R G IR+ such that for any initial condition (So, Qo) G -R, the 
solution converges exponentially to i?o, with a rate rj = 7 A y. We now wish to obtain a 
concentration result for the perturbed system (fT2"j) . that is give a proof of Theorem 11.31 
To this aim, we shall divide our proof in several subsections. 

Notation 3.1. We will set Z\ for the couple (S^,Q e t ), and Z® for the solution to the 
deterministic equation (TJ2J). 

3.1. Reduction of the problem. Recall that Theorem 1 1 . 3 1 st at es an exponential bound 
(valid for p small enough) of the form 

P (\\Z £ - Eo\\oo,i > 2p) < exp f-^J^) , (13) 

on any interval of the form I = [k\ ln(c/p)/r/; K 2 ln(c/p)/?7] and 1 < k\ < k 2 < K3 such 
that A > K3/77. 
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A first step in this direction is to consider a generic interval of the form / = [a, b], and 
write 

P (\\Z £ - EoW^j > 2p) = P - £ |L,f > 2p) n (\\Z° - EoW^j > p)) 

+ p [(\\z £ - ^ |L,f > 2 P ) n (||z° - EoW^j < P )) , 

which yields 

p(||Z £ -£ |L,f>2p) < A 1 + A 2l 

with 

A 1 = p(\\Z°-E \\.>p), and A 2 = p(\\Z £ -Z°\\ f >p). (14) 



Moreover, the term A\ is easily handled: owing to fjlip . we have A\ = as soon as 
a = Ki ln(c/p)/?7 with K\ > 1. In order to prove f|T3|) . it is thus sufficient to check the 
following identity: 



\Z* - Z ^; > p) <exp f-^J^) 



(15) 



on any interval of the form / = ln(c/p)/?7; ln(c/p)/ry] and 1 < < k 2 < We 
shall focus on this inequality in the next subsection. 

3.2. Exponential concentration of the stochastic equation. We will now give a 
general concentration result for Z £ — Z° on suitable time scales as follows: 

Proposition 3.2. Let Z e be the solution to [W\l . Then there exists Eq = Eq(M,t) such 
that, for any p < 1 and e < e we have 



|^-Z°m [0 , T] >p)<exp ("^) 



(16) 



where c 2 , k 2 are strictly positive constants which do not depend on p,E, but both depend 
on our set of parameters a, k, a, d, m, b, M. 

Proof. For notational sake, let us abbreviate ||/||oo,[o,r] into ||/||oo throughout the proof. 
In order to bound Z e — Z°, we first seek a bound for S £ — S°. To this aim we notice 
that for the deterministic function S° and thanks to relation fllip . one can find a constant 
K i = Ki(a,k,a,d,m,b) such that US' Hoc < K\. Set also J\ := J cr(S £ ) odW}. Then 



1£ OO 



|Sf-S, u | < 



< 



/ I (a - ka(QD) S £ s -{a- ka(Q° s )) S°| ds + e \j}\ 
Jo 







ds + 


/ k \ 




Jo 



< [ (a + kM)\S £ s - S° s \ds + Kl k [ \Q £ s -Q° s \ds + E\Jl\. (17) 
Jo Jo 

Analogously, setting J 4 2 := J^a(Q £ s ) o dW^, we obtain 

\Qt-Qt\< [ (m + k(b-l)K 1 )\Q £ s -Q° s \ds + [ k{b- 1)M\S £ S - S° s \ds + e\J?\. (18) 
Jo Jo 
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Hence, putting together (|17|) and (ITS]) , we get the existence of two positive constants 
K2, k 3 such that 

-t 

" Ms, 



|zf-z t T<« 2 s 2 (|j t 1 | 2 + |j t 2 | 2 ) + « 3 / l^ £ -^. 

Jo 



012, 



and by a standard application of Gronwall's lemma, we get for all t G [0, r]: 

|^-Z°| 2 < ^[iJ/P + l^l^expM) 

< K2e 2 [\J}\ 2 + \J?\ 2 ]exp(KsT). (19) 
Let us now go back to our claim (TIB"]) : thanks to inequality (|19]l . we have 

p(l|z=-z»|U> P )=p(||z"-z»||L>p 2 ) 



with 



T 1 = P[\\J 1 \\ 00 > K f - ) , and T 2 = P(||J 2 || 0O >^^ 

V eexp(K 5 r)y V eexp(K 5 r)y 

We now proceed to bound the quantity Ti, and to this aim we first write J\ in terms 
of Ito's integrals: according to [BJ Definition 3.13 p. 156], 



where (•, •) stands for the bracket of two semi-martingales. Invoking equation ( TT2j) and 
ordinary rules of Stratonovich differential calculus, it is also readily checked that 

ft 



Jo 



where V is a process with bounded variation. We thus end up with the expression J\ 
Ml + V t l , where 

M l t = j\{St)dWl and V t l = | jf aa'(S!)ds, 
and decompose T x accordingly into T\ < T^i + T 12 , with 

r ltl = pf||M 1 || eo > 7 . V and r 1 , 2 = pf||v 1 || 0O > "f 

V eexp(K 3 r)/ V eexp(K 3 r) 

We now bound the terms T^i and Ti. 2 separately. 

The term T 12 is easily bounded thanks to some deterministic arguments. Indeed, since 
aa'(x) < C(M + 1) for any x G R+, we have H^Hoo < C(M + l)er, so that for any 
p < 1 and e < £i := (k 4 /(C(M + l)r exp(K 3 r))) 1 / 2 , we have T 12 = 0. As far as T^i is 
concerned, one can apply the exponential martingale inequality (see, for instance, |3|) for 
stochastic integrals in order to get 

/ K.r 2 

Ti i < exp 



K 4 p 



M 2 exp(K S r)s 2 
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Putting together the estimates for and T\ )2 , we have thus obtained 

T\ < exp -— — , 

\ M 2 ex.p(K 3 T)e 2 J 

for any p < 1 and e < S\ := (k a /(C(M + l)r exp(K 3 r))) 1 / 2 . We let the reader check that 
the term T 2 can be handled along the same lines, which finishes our proof. 

□ 



3.3. Deviation from equilibrium. Let us now prove inequality (1131) : recall that we 
have decomposed P(||Z e — -Eolloo,/ > 2p) into A\ + A 2 defined by fU4p . Furthermore, 
A\ = when I is of the form [a, b] with a = K\ ln(c/p)/r). 

In order to complete our result, let us analyze the term A 2 in the light of inequality ffl"6j) . 
Indeed, in order to go from f lT6|) to f|T5|) . it is sufficient to choose p, r, A such that 

p 2 exp(-K 2 r) > p 2+A , 

which is achieved for r < 6 := X\n(l/p)/K2- Hence our claim is satisfied on the interval 
/ = [a, b]. We now have to verify that this interval is nonempty, namely that a < b. This 
gives a linear equation in ln(l/p), of the form 

^ [ln(l/p) + ln(c)] < -ln(l/p). 

and the reader might easily check that the following conditions are sufficient: 
(i) The linear terms satisfy — < — , that is A > — 



(ii) We take p small enough, namely p < po in order to compensate the term ln(c). 
The proof of ( !T3|) is now finished. 

3.4. Extension to the delayed system. Let us deal now with the delayed case: as 
mentioned in the introduction, we consider the system 

dSt = [a- ka(Qt)} Stdt + ea(St) o dW/ 

dQl =[d- mQl - ka(Ql)St + kbe~^a(Ql_ c )Sl ( ] dt + ea(Q £ t ) o dWl 

where for any t e [-(, 0] and for any e > 0, (S%, Qf) = {S?, Qf). 

Under Hypothesis II. 1[ [2~Tl and l2~7l we have shown the existence of a unique equilibrium 
i?o for the deterministic system (pQ), corresponding to (l20|) with e = 0. Following the 
non-delayed case, we wish to obtain a concentration result for the perturbed system (j20p . 
as is given in Theorem II. 51 

The proof of this result can be carried out almost exactly as for Theorem 11.31 Let us 
only point out the main difference: how to get an equivalent of inequalities ()17|) and (j!8|) . 
To this aim, we set again J\ := J Q cr(Sg) o dW} and J t 2 := J* Q * ea(Q e s ) o dW^. Then in the 
delayed case, relations (|17|) and (|18|) become 

\S £ t - Sf\ < [ {a + kM)\S e s - S° s \ds + Kl k [ \Q £ S - Q° s \ds + e\Jl\, (21) 
Jo Jo 
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and 



\Qt-Qt\ < I ( m + kKi)\Q e s - Q° s \ds + f kM\S £ s - S° s \ds + e\J?\ (22) 
Jo Jo 

+ f kbMe~^\S £ s _ c - S° 3 _ c \ds+ [ kbhe-^\Q e s _ c -Q° s _ c \ds. 
Jo Jo 

Using that for any t G [— C, 0] and for any e > 0, (Sf,Q e t ) = (S®,Qt) we can write the 
bounds 

rt-t r t 
kbMe-i* \S £ S _ C - S° s _ ( \ds = / kbMe~ K \S e s - S° s \ds < / kbMe~^\S £ s - S° s \ds, 

Jo Jo 
rt-C ft 
kbk ie -^\Q £ s _ ( - Q° s _ ( \ds = / kbhe-^lQl - Q° s \ds < / kbhe^lQl - Q° s \ds 

Jo Jo 

Then, putting these last bounds in (1211) and (1221) we get the existence of two positive 
constants /t 2 , k 3 such that 

\z° - z°| 2 < k 2 {\j l t \ 2 + |j t 2 | 2 ) + k 3 f \zt - Z^ds. 

Jo 

From this point, the proof follows exactly as for Theorem 11.31 



J o 



4. Numerical simulations 

This final section is devoted to a presentation of some numerical simulations for the sys- 
tem described by equation ([2]). We have chosen the parameters (a, k, d, m, b, according 
to some real data observed by the Molecular Biology Group of the Department of Genetics 
and Microbiology at Universitat Autdnoma de Barcelona. However, since the stochastic 
effects we are dealing with cannot be taken into account in laboratory experiments, we 
did not try to match existing curves like those of [T]. We have chosen instead to compare 
theoretical and noisy dynamics in order to see that the quantities S and Q are close to 
their equilibrium after a reasonable amount of time, in spite of randomness. 

It is worth noticing at this point that the parameters we have chosen for our simulations 
do not meet the conditions stated at Hypothesis 12.71 Indeed, those conditions were 
imposed in order to obtain our theoretical large deviations type results with a reasonable 
amount of efforts, but might be too restrictive to fit to real data experiments. Nevertheless, 
our simulations turn out to be satisfactory, since we observe that the solution (St, Qt) 
converges to Eq for small values of e in a reasonable amount of time, regardless of the 
violation of Hypothesis 12.71 

Specifically, we have simulated trajectories with parameters estimated on an experiment 
involving Salmonella ATCC14028 bacteria and UAB_Phi78 virus. From the experiments 
conducted by the mentioned group we have chosen the parameters as: 

(a, k, d, m, b, () = (12.1622, 27.36, 0.1, 0.1947, 61, 0.01875). 

We have also put M — 10, fi — 0.5, and we have taken the initial conditions Sqj = 
4.8e a< -* + ^, Qo t t = for t E [C?0]. The time is expressed in days and the amount of virus 
and bacteria are expressed in tens of millions of units. 
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Figure 1. Simulation of the trajectories of S and Q with real parame- 
ters for the Salmonella ATCC14028 bacteria and UAB_Phi78 virus for the 
deterministic case (e = 0), for e = —3 (red curve) and e — 1 (blue curve). 

Our simulations are summarized at Figure [U in which different paths of the processes S 
and Q are computed. We have used an Euler type discretization scheme for our equations, 
implemented with the R software. We have then plotted the deterministic case (e = 0) 
plus the curves corresponding to several values of e (namely e = —3, 1). As mentioned 
before, the fluctuations of S and Q (which are obviously due to the randomness we have 
introduced) do not prevent them to converge to equilibrium. 
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